! Copyright (c) 2021~2023, ZUO Zhihua.
!> random number
module random_module

#ifdef REAL64
    use, intrinsic :: iso_fortran_env, only: rk => real64
#else
    use, intrinsic :: iso_fortran_env, only: rk => real32
#endif
    implicit none

    private
    public :: randn, randu, rk

    real(rk), parameter :: pi = acos(-1.0_rk)

contains

    !> uniform distribution
    impure elemental subroutine randu(x, min, max)
        real(rk), intent(in) :: min, max
        real(rk), intent(out) :: x

        call random_number(x)
        x = min + x*(max - min)

    end subroutine randu

    !> normal distribution, Box-Muller method
    impure elemental subroutine randn(x, mean, std)
        real(rk), intent(in) :: mean, std
        real(rk), intent(out) :: x
        real(rk) :: y

        call random_number(y)
        call random_number(x)
        x = mean + std*sqrt(-2.0_rk*log(y))*cos(2.0_rk*pi*x)

    end subroutine randn

end module random_module
